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Abstract 


Observational  data  from  the  Global  Temperature  and  Salinity  Profile  Program  (GTSPP) 
were  used  to  calculate  the  upper  ocean  heat  content  (OHC)  anomaly.  The  thickness  of  the  upper 
layer  is  taken  as  300  m  for  the  Pacific/ Atlantic  Ocean  and  150  m  for  the  Indian  Ocean  since  the 
Indian  Ocean  has  shallower  thermoclines.  First,  the  optimal  spectral  decomposition  (OSD) 
scheme  was  used  to  build  up  monthly  synoptic  temperature  and  salinity  data  set  for  January  1990 
to  December  2009  on  l°xl°  grids  and  the  same  33  vertical  levels  as  the  World  Ocean  Atlas. 
Then,  the  monthly  varying  upper  layer  OHC  field  (H)  was  obtained.  Second,  a  composite 

analysis  was  conducted  to  obtain  the  total-time  mean  OHC  field  ( H )  and  the  monthly  mean 
OHC  variability  ( H  ),  which  is  found  an  order  of  magnitude  smaller  than  H  .  Third,  an  empirical 

orthogonal  function  (EOF)  method  is  conducted  on  the  residue  data  (H  ),  deviating  from  H  +  H  , 
in  order  to  obtain  interannual  variations  of  the  OHC  fields  for  the  three  oceans.  In  the  Pacific 
Ocean,  the  first  two  EOF  modes  account  for  51.46%  and  13.71%  of  the  variance,  representing 
canonical  El  Nino/La  Nina  (EOF-1)  and  pseudo  El  Nino/La  Nina  (i.e.,  El  Nino  Modoki)  (EOF-2) 
events.  In  the  Indian  Ocean,  the  first  two  EOF  modes  account  for  24.27%  and  20.94%  of  the 
variance,  representing  basin-scale  cooling/wanning  (EOF-1)  and  Indian  Ocean  Dipole  (EOF-2) 
events.  In  the  Atlantic  Ocean,  the  first  EOF  mode  accounts  for  49.26%  of  the  variance, 
representing  a  basin-scale  cooling/warming  (EOF-1)  event.  The  second  EOF  mode  accounts  for 
8.83%  of  the  variance.  Different  from  the  Pacific  and  Indian  Oceans,  there  is  no  zonal  dipole 
mode  in  the  tropical  Atlantic  Ocean.  Fourth,  evident  lag  correlation  coefficients  are  found 
between  the  first  principal  component  of  the  Pacific  Ocean  and  the  Southern  Oscillation  Index 
with  a  maximum  correlation  coefficient  (0.68)  at  1 -month  lead  of  the  EOF-1  and  between  the 
second  principal  component  of  the  Indian  Ocean  and  the  Dipole  Mode  Index  with  maximum 


2 


values  (around  0.53)  at  1-2  month  advance  of  the  EOF-2.  It  implies  that  OHC  anomaly  contains 
climate  variability  signals. 


Keywords:  Dipole  Mode  Index  (DMI),  El  Nino;  pseudo-El  Nino;  Global  Temperature  and 
Salinity  Profile  Program  (GTSPP);  Indian  Ocean  dipole;  monthly  varying  global  (T,  S)  field; 
optimal  spectral  decomposition  (OSD);  Southern  Oscillation  Index  (SOI);  dipole  mode  index 
(DMI);  upper  ocean  heat  content 
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1.  Introduction 


In  the  past  decade,  several  new  coupled  ocean-atmosphere  phenomena,  such  as  the 


pseudo-El  Nino  (or  sometimes  called  El  Nino  Modoki)  and  the  Indian  Ocean  Dipole  (IOD),  were 


discovered  and  recognized  important  in  climate  variability.  The  pseudo-El  Nino  is  characterized 
by  warmer  sea  surface  temperature  anomalies  (SSTA)  in  the  central  equatorial  Pacific  and  cooler 
SSTA  in  the  eastern  and  western  equatorial  Pacific  (Ashok  et  al.  2007;  Weng  et  al.  2007),  which 
is  different  from  the  El  Nino  with  anomalous  warming  in  eastern  equatorial  Pacific.  El  Nino  and 
pseudo-El  Nino  have  different  teleconnection  patterns.  Taking  the  Atlantic  Ocean  as  an  example, 
less  tropical  storms  and  hurricanes  occur  during  El  Nino  events  (see  the  website: 
http ://ww2 010. atmos . uiuc . edu/ ( Gh)/ guides/ mtr/hurr/ enso . rxml) ;  and  more  tropical  storms  and 
hurricanes  appear  during  pseudo-El  Nino  events  (Weng  et  al.  2007).  Ashok  et  al.  (2007)  and 
Weng  et  al.  (2007)  described  in  detail  the  difference  in  teleconnections  and  climate  variability 
between  the  pseudo-El  Nino  event  in  2004  and  El  Nino  events.  The  IOD  was  identified  by 
Saji  et  al.  (1999)  and  Vinayachandran  et  al.  (1999). 

Climate  systems  are  complex.  One  way  to  simplify  the  climate  systems  is  to  represent 
climate  variability  of  atmospheric  circulations  by  teleconnection  patterns,  such  as  the  Southern 
Oscillation  Index  (SOI),  representing  the  El  Nino/La  Nina  events,  and  the  Dipole  Mode  Index 
(DMI),  representing  the  IOD  events.  The  SOI  was  first  to  utilize  equivalent  barotropic  seesaw  in 
atmospheric  pressure  between  the  southeastern  tropical  Pacific  and  the  Australian-Indonesian 
regions  (Walker  and  Bliss  1937).  A  popular  formula  for  calculating  the  monthly  SOI  is  proposed 
by  the  Australian  Bureau  of  Meteorology: 


diff 


(1) 
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Here, /?diff  is  the  mean  sea  level  pressure  of  Tahiti  minus  that  of  Dawin  for  that  month;  pdiff  'j  is 

the  long-term  average  of  pm  for  the  month  in  question;  and  SD(pdiff)  is  the  long-term  standard 
deviation  of  pm  for  the  month  in  question.  Positive  values  of  the  SOI  (Fig.  la)  are  associated 
with  stronger  Pacific  trade  winds  and  warmer  sea  temperatures  to  the  north  of  Australia. 

The  dipole  mode  index  (DMI)  (Fig.  lb)  was  introduced  by  Saji  et  al.  (1999)  to  define  the 
strength  of  the  IOD  events.  It  is  defined  by  the  difference  in  SSTA  between  the  tropical  western 
Indian  Ocean  (10°S-10°N,  50°-70°E)  and  the  tropical  southeastern  Indian  Ocean  (10°S-equator, 
90°-110°E).  During  positive  IOD  events  (i.e.,  positive  DMI),  a  wann  SSTA  occurs  in  the 
western  equatorial  Indian  Ocean  and  a  cold  SSTA  appears  in  the  eastern  equatorial  Indian  Ocean 
(“La  Nina-type”  in  the  Indian  Ocean),  which  causes  an  easterly  wind  anomaly,  and  in  turn 
weakens  the  eastward  equatorial  jet  (EEJ).  During  negative  DMI  (i.e.,  negative  IOD  events),  a 
cold  SSTA  occurs  in  the  western  equatorial  Indian  Ocean  and  a  warm  SSTA  appears  in  the 
eastern  equatorial  Indian  Ocean  (“El  Nino-type”  in  the  Indian  Ocean),  which  causes  a  westerly 
wind  anomaly,  and  in  turn  enhances  the  eastward  equatorial  jet.  Chu  (2010)  confirmed  the 
relationship  between  IOD  and  eastward  equatorial  jet  using  observational  data. 

Upper  ocean  heat  content  (OHC)  has  been  recognized  important  for  climate  variability. 
For  example,  Hasegawa  and  Hanawa  (2003)  shows  the  existence  of  anticlockwise  propagation 
of  OHC  anomalies  in  the  tropical  Northern  Pacific,  which  is  related  closely  to  the  ENSO  events. 
However,  the  pseudo-El  Nino  and  IOD  events  were  identified  as  ocean-atmospheric  coupled 
system  from  ocean-atmospheric  data  analysis  with  SST  rather  than  OHC  as  a  key  variable  for  the 
ocean.  Questions  arise:  Can  the  three  major  tropical  systems  (ENSO,  pseudo-El  Nino  and  IOD) 
be  identified  from  the  upper  layer  OHC  field  alone?  What  are  connections  between  OHC 
anomaly  and  climate  indices?  To  answer  these  questions,  an  observational  dataset  from  the 
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Global  Temperature  and  Salinity  Profile  Program  (GTSPP)  was  used.  GTSPP  has  been 
developed  since  1990  with  international  cooperation  (see  the  website:  http://www.iode.org/). 

With  more  observed  ( T ,  S)  profile  data  especially  with  the  Argo  floats,  it  is  possible  to 
establish  gridded  (T,  S)  dataset  with  a  sufficient  resolution  in  space  and,  especially,  in  time  using 
the  optimal  spectral  decomposition  (OSD)  method  (Chu  et  al.  2003a,  b).  From  the  gridded 
dataset,  it  is  able  to  compute  the  upper  ocean  heat  content  (OTIC).  The  thickness  of  the  upper 
layer  is  taken  as  300  m  for  the  Pacific/ Atlantic  Ocean  and  150  m  for  the  Indian  Ocean  since  the 
Indian  Ocean  has  shallower  thermoclines.  Composite  were  used  to  get  gridded  data  of  OHC 
seasonal  variability  and  anomaly. 

This  paper  presents  such  an  effort  and  is  outlined  as  follows.  Section  2  depicts  the 
GTSPP  dataset,  Section  3  presents  characteristics  of  spatially  (1°  resolution)  and  temporally  (1 
month  increment)  varying  OHC  (H:  H30o  for  the  Pacific/ Atlantic  Oceans  and  Hi50  for  the  Indian 
Ocean)  after  using  the  OSD  method,  a  composite  analysis  to  get  the  OHC  seasonal  variability, 
and  an  EOF  analysis  to  obtain  the  OHC  anomaly.  Sections  4-6  describe  temporal  and  spatial 
variability  of  OHC  anomaly  and  the  connection  to  the  climate  variability.  Section  7  lists  the 
conclusions. 

2.  GTSPP 

GTSPP  is  a  cooperative  international  project.  It  seeks  to  develop  and  maintain  a  global 
ocean  temperature-salinity  resource  with  data  that  are  both  up-to-date  and  of  the  highest  quality 
possible.  Making  global  measurements  of  ocean  temperature  (7)  and  salinity  ( S)  quickly  and 
easily  accessible  to  users  is  the  primary  goal  of  the  GTSPP.  Both  real-time  data  transmitted  over 
the  Global  Telecommunications  System,  and  delayed-mode  data  received  by  the  National 
Oceanographic  Data  Center  (NODC)  are  acquired  and  incorporated  into  a  continuously  managed 
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database.  The  GTSPP  handles  all  temperature  and  salinity  profile  data.  This  includes 
observations  collected  using  water  samplers,  continuous  profding  instruments  such  as  Argo 
floats,  CTDs,  thermistor  chain  data  and  observations  acquired  using  thermosalinographs.  The 
GTSPP  went  through  quality  control  procedures  that  make  extensive  use  of  flags  to  indicate  data 
quality.  To  make  full  use  of  this  effort,  participants  of  the  GTSPP  have  agreed  that  data  access 
based  on  quality  flags  will  be  available  (Sun  et  al.  2009).  That  is,  GTSPP  participants  will  permit 
the  selection  of  data  from  their  archives  based  on  quality  flags  as  well  as  other  criteria.  These 
flags  are  always  included  with  any  data  transfers  that  take  place.  Because  the  flags  are  always 
included,  and  because  of  the  policy  regarding  changes  to  data,  as  described  later,  a  user  can 
expect  the  participants  to  disseminate  data  at  any  stage  of  processing.  Furthermore,  GTSPP 
participants  have  agreed  to  retain  copies  of  the  data  as  originally  received  and  to  make  these 
available  to  the  user  if  requested  (Sun  et  al.  2009).  Interested  readers  are  referred  to  the  following 
website  (http://www.iode.org/)  of  the  International  Oceanographic  Commission  of  UNESCO. 

There  are  total  6,662,209  GTSPP  stations  with  temporally  increasing  data  (Fig.  2).  The 
yearly  number  of  ( T ,  S)  profdes  is  less  than  80,000  from  1990  to  1999.  Since  2000,  as  the  Argo 
floats  into  practice,  the  yearly  number  of  ( T ,  S)  profiles  increases  exponentially  close  to  1.7 
million  in  2009.  The  GTSPP  stations  are  over  the  global  ocean  (Fig.  3).  The  OSD  method  (see 
Appendix  A)  is  used  to  establish  3D  global  gridded  ( T ,  S)  dataset  with  monthly  increment,  l°xl° 
horizontal  resolution,  and  33  vertical  layers  as  the  World  Ocean  Atlas 
(http://www.nodc.noaa.gov/OC5/WOAQ5/pr  woa05.html),  which  is  called  as  the  GTSPP 
gridded  data. 

3.  Spatial  and  Temporal  Variability  of  OHC 

3.1.  Gridded  H30o  and  HiS0  Datasets 
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From  the  GTSPP  gridded  dataset,  the  upper  OHC  ( H)  with  the  unit  of  J/m“  is  calculated 
by 

0 

H(r;.,  rk,tl)  =  cp  J  pTdz,  h  =  150  m,  300  m  ,  (2) 

-h 

where  cp  is  the  specific  heat  of  the  seawater  [=  3850  J/(kg  K)];  and  p  is  the  density,  which  is 
calculated  from  the  ( T ,  S)  data  using  the  equation  of  the  state;  the  vectors  r,  represent  locations 
of  2D  horizontal  grids;  i  =  1,  2,  ...,  I,  where  /  is  the  total  number  of  the  horizontal  grid  points; 
ik  =  1990,  1991,  ...,  2009  is  the  time  sequence  in  years;  and  ti  =  1,  2,  ...,  12  is  the  monthly 
sequence  within  a  year. 

3.2.  Total-Time  Mean 

Before  investigating  the  annual  variation  of  OHC,  we  define  the  following  two  temporal 
averages: 

H(r  ,t,)  =  —  VH(r,r,  ,tX  Ar  =  20  (1990  to  2009) ,  (3a) 

Ary 

which  are  the  annual  mean  values,  and 

fiw=-rrZfl<w,).  (3b) 

iz  /= i 

which  are  the  total-time  mean  values.  The  total-time  mean  (1990-2009)  of  global  upper  ocean 
heat  content  with  the  upper  ocean  depth  of  300  m  H,(l0 (r)  (Fig.  4a)  shows  following  features: 
(1)  it  is  usually  less  than  1.2  xlO10  J/m 2  north  of  40°N  and  south  of  40°S.  (2)  The  high 

H300(r/)  areas  are  associated  with  the  subtropical  gyres  in  the  Atlantic  and  Pacific  Oceans.  For 
example,  in  the  Atlantic  Ocean,  a  strong  warm  core  (3.0x10  J/m-)  is  located  in  the  Gulf  Stream 
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region.  It  extends  northeastward  and  reaches  the  northwest  African  Coast.  Another  warm  area 
(2.8><10  J/m  )  is  located  in  the  South  Atlantic  with  a  triangular  shape  from  the  Brazilian- 
Argentine  Coast  (15°S  to  40°S)  to  the  southwest  African  Coast  (around  40°S).  Between  the  two 
warm  regions,  there  is  a  cool  area  (1.9x10  J/m")  including  the  equatorial  region.  In  the  Pacific 
Ocean,  a  warm  core  in  the  western  North  Pacific  (the  Kuroshio  and  its  extension  regions)  and 
the  other  in  the  western  equatorial/South  Pacific  with  the  values  around  3.0X10  J/m  .  The 
eastern  Pacific  has  less  values  around  2.0X10  J/m-.  The  total-time  mean  of  global  upper  ocean 

heat  content  with  the  upper  ocean  depth  of  150  m  H150(r(.)  (Fig.  4b)  shows  similar  features  with 
reduced  values  of  H(r(.) . 

3.3.  Mean  Seasonal  Variability 

The  annual  mean  values  deviated  from  the  total-time  mean,  H ( r) , 

H(r„()  =  H(ri,()-fi(r,),  (4) 

lead  to  composite  features  of  the  mean  seasonal  variability.  The  mean  seasonal  variability  of 
OHC  with  upper  ocean  thickness  of  300  m  (Fig.  5a)  and  150  m  (Fig  5b)  shows  the  following 
features:  (l)H300(r(,() ,  with  maximum  values  around  2.6x10  J/m  ,  is  generally  an  order  of 

magnitude  smaller  than  the  total-time  mean  H300(r(.) .  Such  a  feature  also  exists  in  Hl50(r,() . 

(2)  Its  positive  values  show  eddy-like  structure  and  negative  values  are  more  spatially  uniform. 

(3)  Usually,  negative  (positive)  values  occur  in  the  northern  (southern)  hemisphere  in  boreal 
winter  (upper  panels  in  Fig.  5),  and  positive  (negative)  values  occur  in  the  northern  (southern) 
hemisphere  in  boreal  summer  (lower  panels  in  Fig.  5). 

3.4.  OHC  Anomaly 
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OHC  anomaly  is  the  deviation  from  its  annual  mean  value 


H(r,.  ,rk,t,)  =  H(r,  ,rk,t,)-  H(r ,  t, ) ,  (5) 

which  is  re-arranged  into  a  NxP  matrix,  H(r ),t  ),  p  =  1,  2,  . P.  Here  P  =  240  is  the  total 

number  of  months  (i.e.,  total  time  points).  The  empirical  orthogonal  function  (EOF)  analysis 
(Lorenz  1956)  was  conducted  on  the  heat  content  anomaly  data 

H  =  ((,)£, (r,),  (6) 

j 

where  E}(r\)  is  the  jth  EOF  mode  with  the  unit  of  J/m  ;  and  PCj{tp)  is  the  jth  principal 
component  without  unit.  Here,  j  =  1,  2,  . . .,  J;  Jis  truncated  number  of  the  EOF  modes.  The  EOF 
analysis  was  conducted  for  the  individual  ocean  such  as  on  H300(r.,?  )  for  the  Pacific  and 

Atlantic  Oceans  and  on  H150(r.,f  )  for  the  Indian  Ocean. 

4.  Pacific  Ocean 

Table- 1  shows  the  variance  of  H300(r.,f  )  for  the  first  six  EOFs  with  EOF- 1  containing 

51.46%  of  variance  and  EOF-2  containing  13.71%  of  variance.  Here,  the  effort  is  concentrated 
on  analyzing  the  first  2  EOF  modes. 

4.1.  EOF-1  Mode 

The  EOF-1  mode  (Fig.  6)  shows  the  positive  phase  of  a  zonal  dipole  structure  (i.e., 
positive-negative):  The  western  Pacific  (near  half  of  the  Pacific  Ocean)  is  positive  with  a  strong 
maximum  center  located  at  the  western  tropical  Pacific  warm  pool  region  north  of  New  Guinea 
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Q  2, 

with  maximum  value  higher  than  2.0x10  J/rrr.  The  eastern  Pacific  is  negative  with  minimum 
value  lower  than  -0.8  xlO9  J/m2  located  in  the  east  equatorial  Pacific. 

The  first  principal  component  for  1 990-2009, PCx(tp),  fluctuates  between  -1  to  1,  which 

makes  the  magnitude  of  the  EOF-1  pattern  around  2.0x10  J/m",  comparable  to  the  magnitude  of 
the  seasonal  variability  (2.6xl09  J/m2).  It  has  an  evident  upward  trend  from  PCX (t p)  generally 

negative  before  1998  to  PCx(tp)  generally  positive  after  1999  except  2002  (Fig.  7a). 
ForPC,^)  >  0,  the  horizontal  variability  of  OHC  anomaly  shows  the  positive  phase  (Fig.  6a), 
representing  the  Fa  Nina  events.  ForEC^)  <  0 ,  the  horizontal  variability  of  OHC  anomaly 
shows  the  negative  phase,  representing  the  El  Nino  events.  The  values  of  PCx(tp)  are  generally 
negative  in  1990-1998.  During  this  period,  PCx(tp)  has  minimum  value  in  1991-1992,  1993, 

1995,  and  1997-1998.  They  correspond  to  1991-1995  long  El  Nino  event,  and  1997-1998  El 
Nino  event.  The  values  of  PCx(tp)  are  positive  in  1999-2009,  which  implies  that  there  is  no  El 

Nino  event  (in  the  classical  sense)  during  this  period.  Taking  1997-1998  El  Nino  event  as  an 
example,  PCx(tp)  is  -1.0  in  November  1997.  The  OHC  has  wann  anomaly  greater  than  0.80x  109 

J/nr  in  the  eastern  Pacific,  and  cold  anomaly  in  the  western  Pacific. 

4.2.  EOF-2  Mode 

The  EOF-2  mode  (Fig.  6b)  shows  the  positive  phase  of  a  different  zonal  structure:  Small 
portion  of  the  western  Pacific  Ocean  is  positive  with  a  maximum  (around  0.45x10  J/m")  located 
near  the  Indonesian  islands.  The  rest  of  the  Pacific  Ocean  is  negative  with  a  minimum  value 
around  -l.OxlO9  J/m 2  appearing  in  the  equatorial  central  Pacific  (130°W  -  120°W,  3°N  -  3°S). 
The  second  principal  component,  PC2{tp) ,  fluctuates  between  -1  to  1  during  1990-2009, 
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which  makes  the  magnitude  of  the  EOF-2  pattern  around  1.0x10  J/m",  a  little  smaller  than  the 
magnitude  of  the  seasonal  variability  (2.6x10  J/m").  Different  from  PC](tp) ,  the  second 

principal  component,  PC2 (t  )  ,  does  not  have  a  trend  (Fig.  7b).  For PC2(t  )  >  0  ,  the  horizontal 

variability  of  OHC  anomaly  shows  the  positive  phase,  representing  the  Fa  Nina  event. 
For  PC2(tp)  <  0,  the  horizontal  variability  of  OHC  anomaly  shows  the  negative  phase  (i.e., 

opposite  as  shown  in  Fig.  6b).  That  is,  negative  OHC  anomalies  occupy  the  small  portion  of  the 
western  Pacific  Ocean;  positive  OHC  anomalies  occupy  rest  of  the  oceans  with  warmest 
anomaly  in  the  central  Pacific  Ocean,  representing  the  pseudo-El  Nino  events.  The  values  of 
PC2(tp)  were  negative  in  short  periods  in  1992  ,  1993,  and  1996,  January  1997  to  July  1998, 

March  2002  to  February  2007,  and  the  whole  year  of  2009.  Among  them,  the  period  of  January 
1997  to  July  1998  is  very  special,  when  both  PCx(tp)  and  PC2(t  )  are  negative,  which  implies 

co-existence  of  an  El  Nino  event  from  the  EOF-1  and  a  pseudo-El  Nino  event  from  the  EOF-2. 
Furthermore,  the  value  of  PC2(tp)  is  -1.0  in  October  1997.  The  OHC  of  the  EOF-2  mode  has 

warm  anomaly  greater  than  1.0X109  J/m2  in  the  central  Pacific;  cold  anomalies  in  the  western 
Pacific  with  a  minimum  around  -0.45 x  10  J/m"  near  Indonesian  islands. 

4.3.  Phase  Space  Trajectory 

Temporal  variation  of  the  two  modes  can  also  be  represented  by  the  phase  diagram  with 
PC]  as  the  horizontal  axis  and  PCi  as  the  vertical  axis  (Fig.  8).  The  trajectory  in  the  phase  space 
starts  from  January  1990  to  December  2009  (a  star  enclosed  by  a  circle).  The  number  enclosed 
by  a  circle  denotes  January  of  the  year  (90  means  1990,  ...);  the  arrows  show  the  time 
advancement;  the  red  color  refers  to  negative  values  of  SOI  and  blue  color  refers  to  positive 
values  of  SOI.  The  trajectory  represents  evolution  of  the  OHC  anomaly  in  various  phases 
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relating  to  climate  variability.  It  is  noted  that  both  positive  (blue)  and  negative  (red)  SOI  are 
found  in  all  four  quadrants. 

For  the  trajectories  between  two  consecutive  years  in  the  quadrant  of  (+,  +)  such  as  99- 
00-01  (representing  January  1999  to  January  2001),  08-09  (representing  January  2008  to  January 
2009),  both  EOF-1  and  EOF-2  modes  are  in  the  positive  phase,  implying  strong  La  Nina  in  the 
tropical  Pacific.  For  the  trajectories  between  two  consecutive  years  in  the  quadrant  of  (-,  -)  such 
as  March  97-  January  98,  both  EOF-1  and  EOF-2  modes  are  in  the  negative  phase,  which  implies 
co-existence  of  El  Nino  and  pseudo-El  Nino  in  the  tropical  Pacific.  For  the  trajectories  between 
two  consecutive  years  in  the  quadrant  of  (+,  -)  such  as  04-05-06-07,  and  most  months  in  2009, 
the  EOF-1  mode  is  in  the  positive  phase  and  the  EOF-2  mode  is  in  the  negative  phase,  which 
implies  pseudo-El  Nino  in  the  tropical  Pacific.  For  the  trajectories  between  two  consecutive 
years  in  the  quadrant  of  (-,  +)  such  as  most  months  in  1990-1995,  the  EOF-1  mode  is  in  the 
negative  phase  and  the  EOF-2  mode  is  in  the  positive  phase  or  near  zero  such  as  92  and  95, 
which  implies  El  Nino  in  the  tropical  Pacific. 

4.4.  Lag-Correlations  between  (SOI,  DMI)  and  Principal  Components 

Linkage  between  the  OHC  anomaly  in  the  Pacific  Ocean  and  climate  variability  is 
represented  by  the  lag-correlation  coefficients  between  the  (SOI,  DMI)  with  the  principal 
components  PCx{tp)  andPC2(tp) .  The  total  size  of  the  time  series  (n)  is  240  (20  years  of 

monthly  data  between  1990  and  2009)  for  the  no-lag  correlation  coefficients.  For  an  /-month  lag 
correlation,  the  data  size  is  (n  -  f). 

The  lag-correlation  coefficients  show  strong  association  between  PC^t  )  and  s(tp)  with 

a  maximum  values  around  0.68  with  1 -month  advance  of  PC\  (Fig.  9a,  upper  panel),  weak 
association  between  PC2(tp)  and  s(t )  with  maximum  values  around  0.32  with  no-lag  to  3- 
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month  lag  of PC2(t  )  to  s(tp)  (Fig.  9a,  lower  panel),  and  no  association  between  DMI  and 

(PC i,  PC 2)  (Fig.  9b).  This  suggests  that  the  EOF-1  mode  is  strongly  related  to  the  Southern 
Oscillation  with  a  maximum  correlation  at  1 -month  lead  of  the  EOF-1;  and  the  EOF-2  mode 
(pseudo  El  Nino)  is  related  weakly  to  the  Southern  Oscillation  with  a  maximum  correlation  at 
no-lag  to  3 -month  lag  of  the  EOF-2. 

5.  Indian  Ocean 

Table-2  shows  the  variance  of  H150(r.,?  )  for  the  first  six  EOFs  with  EOF-1  containing 

24.27%  of  variance  and  EOF-2  containing  20.94%  of  variance.  Here,  the  effort  is  concentrated 
on  analyzing  the  first  2  EOF  modes. 

5.1.  EOF-1  Mode 

The  EOF-1  mode  (Fig.  10a)  shows  a  basin-scale  cooling/warming  (i.e.,  basin-scale  mode) 
with  the  positive  phase  depicted  as  the  basin-scale  cooling  with  a  strong  cooling  zone  stretching 
from  northeast  of  Madagascar  to  the  west  coast  of  Sumatera  with  minimum  value  lower  than 
-0.45x10  J/m  in  the  west  near  the  Mascarene  Plateau.  The  cooling  weakens  away  from  this 
strong  cooling  zone  with  value  around  -O.lxlO9  J/m2  near  30°S. 

The  first  principal  component  for  1 990-2009, PCx(tp) ,  fluctuates  between  -1  to  1,  which 

makes  the  amplitude  of  the  EOF-1  pattern  around  0.45 x109  J/m 2  (Fig.  11a).  This  value  is  only 
17%  of  the  magnitude  of  the  seasonal  variability  (2.6x10  J/m").  Moreover,  PCx(t)  has  an 

evident  downward  trend  from  generally  positive  values  (basin-scale  cooling)  before  2000  (except 
1997-1998)  with  a  maximum  value  of  +1.0  in  November  1992  to  generally  negative  values 
(basin-scale  warming)  after  2000  with  a  minimum  value  of  -0.98  in  November  2007. 


14 


5.2.  EOF-2  Mode 


The  EOF-2  mode  (Fig.  10b)  shows  the  positive  phase  of  a  zonal  dipole  structure  (i.e., 
positive-negative):  The  western  Indian  Ocean  (near  half  of  the  Indian  Ocean)  is  positive  with  a 
warm  center  (higher  than  0.6X10  J/m")  located  northeast  of  Madagascar.  The  eastern  Indian 
Ocean  is  negative  with  cool  center  (lower  than  -0.8  xlO  J/m  )  located  west  of  Sumatera.  This 
pattern  corresponds  to  the  positive  IOD  event,  and  similar  to  the  pattern  obtained  from  the  ocean- 
atmospheric  coupled  mode  (Rao  et  ah,  2002). 

The  second  principal  component,  PC2  (t  )  ,  fluctuates  between  -1  to  1  during  1990-2009, 

which  makes  the  amplitude  of  the  EOF-2  pattern  around  0.8X109  J/m2,  smaller  than  the 
magnitude  of  the  seasonal  variability  (2.6X10  J/m“).  Different  from  PC](tp) ,  the  second 

principal  component,  PC2(t  )  ,  does  not  have  a  trend  (Fig.  1  lb).  For  PC2(t  )  >  0  ,  the  horizontal 

variability  of  OHC  anomaly  shows  the  positive  phase  (i.e.,  the  positive  IOD  event):  The  western 
Indian  Ocean  is  positive  with  a  warm  center  located  northeast  of  Madagascar;  and  the  eastern 
Indian  Ocean  is  negative  with  a  cool  center  located  west  of  Sumatera.  For  .PC,  (?  )  <  0,  the 

horizontal  variability  of  OHC  anomaly  shows  the  negative  phase  (i.e.,  the  negative  IOD  event): 
The  western  Indian  Ocean  is  negative  with  a  cool  center  located  northeast  of  Madagascar;  and 
the  eastern  Indian  Ocean  is  positive  with  a  warm  center  located  west  of  Sumatera,  representing 
the  negative  IOD  event.  Generally,  PC2(t  )  ,  shows  more  frequent  occurrence  of  positive  values 

before  2000  than  after  2000:  1991,  1992,  1993,  1994-1995,  1997-1998,  2003-2004,  and  2006- 
2008.  Among  them,  the  events  of  1994,  1997,  2003,  and  2006  are  significant.  The  period  of 
January  1997  to  July  1998  is  very  special,  the  positive  IOD  event  [positive  PC2(t  )  for  the 
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Indian  Ocean]  is  co-existence  with  El  Nino  event  [negative  PC,  ( tp )  for  the  Pacific  Ocean]  and 
pseudo-El  Nino  event  [negative  PC2(t  )  for  the  Pacific  Ocean]. 

5.3.  Phase  Space  Trajectory 

Temporal  variation  of  the  two  modes  can  also  be  represented  by  the  phase  diagram  with 
PC]  as  the  horizontal  axis  and  PCi  as  the  vertical  axis  (Fig.  10).  The  trajectory  in  the  phase  space 
starts  from  January  1990  to  December  2009  (a  star  enclosed  by  a  circle).  The  number  enclosed 
by  a  circle  denotes  January  of  the  year  (90  means  1990,  ...);  the  arrows  show  the  time 
advancement;  the  red  color  refers  to  negative  values  of  DMI  and  blue  color  refers  to  positive 
values  of  DMI.  The  trajectory  represents  evolution  of  the  OHC  anomaly  in  various  phases 
relating  to  climate  variability.  It  is  noted  that  both  positive  (blue)  and  negative  (red)  DMI  are 
found  in  all  four  quadrants  and  the  phase  space  trajectory  is  more  chaotic  for  the  Indian  Ocean 
(Fig.  13)  than  for  the  Pacific  Ocean  (Fig.  8).  Such  different  behaviors  in  the  phase  space 
trajectories  need  to  be  further  investigated  using  a  numerical  model,  which  is  beyond  the  scope 
of  this  study. 

5.4.  Lag-Correlations  between  (SOI,  DMI)  and  Principal  Components 

Linkage  between  the  OHC  anomaly  in  the  Indian  Ocean  and  climate  variability  is 
represented  by  the  lag-correlation  coefficients  between  (SOI,  DMI)  and  the  principal  components 
PCx(t  )  and PC2(t  )  .  The  lag-correlation  coefficients  show  no  association  between  PCx(tp)  and 

(SOI,  DMI)  (Fig.  14a,  b,  upper  panels),  evident  negative  association  between  PC2(t  )  and  SOI 
with  a  minimum  value  of  -0.57  for  1 -month  lead  of PC2(tp),  and  evident  positive  association 
between  PC2(t  )  and  DMI  with  a  maximum  value  around  0.52  for  2-month  advance  of PC2(t  ) 
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to  s(ip)  (Fig.  14a,  b,  lower  panels).  This  clearly  shows  the  co-existence  of  the  positive  phase  of 
the  EOF-2  mode,  the  negative  phase  of  SOI  (El  Nino  events),  and  the  positive  IOD  events. 

6.  Atlantic  Ocean 

Table-2  shows  the  variance  of  H300(r.,F  )  for  the  first  six  EOFs  with  EOF- 1  containing 

49.26%  of  variance  and  EOF-2  containing  8.83%  of  variance.  Here,  the  effort  is  concentrated 
on  analyzing  the  first  2  EOF  modes. 

6.1.  EOF-1  Mode 

The  EOF-1  mode  (Fig.  15a)  shows  a  positive  phase  of  the  basin-scale  mode  (showing 
basin-scale  cooling)  with  two  strong  and  two  weak  cooling  zones.  The  two  strong  cooling  zones 
with  minimum  values  lower  than  -0.8X109  J/m2  are  located  at  the  tropics  (15°S  -  15°N)  near 
Brazilian  coast  and  high  latitudes  (50°N  -  65°N)  near  Greenland.  The  two  weak  cooling  zones 
are  located  at  the  mid-latitudes  (15°N-50°N)  with  minimum  values  around  -0. 1  x  ]  09  J/m2  in  the 
western  part  near  the  U.S.  coasts  and  around  -0.4x10  J/m“  in  eastern  parts  near  the  northwest 
African  coast.  The  EOF-1  contours  nearly  follow  the  latitudes  in  the  South  Atlantic  Ocean  with 
reduction  of  variability  southward. 

The  first  principal  component  for  1 990-2009, PCx(tp) ,  fluctuates  between  -1  to  1,  which 

makes  the  amplitude  of  the  EOF-1  pattern  around  0.8X109  J/m2  (Fig.  16a).  This  value  is  smaller 
than  the  magnitude  of  the  seasonal  variability  (~2.5X10  J/m",  see  Fig.  5).  Moreover,  PCx(tp) 

has  an  evident  downward  trend  from  generally  positive  values  (basin-scale  cooling)  before  2000 
(except  1997-1998)  with  a  maximum  value  of  +1.0  in  April  1990  to  generally  negative  values 
(basin-scale  warming)  after  2000  except  2002  with  a  minimum  value  of  -0.83  in  March  2005. 
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6.2.  EOF-2  Mode 


The  EOF-2  mode  (Fig.  15b)  shows  the  positive  phase  of  latitudinal  alternative 
cooling/warming  zones:  The  South  Atlantic  is  positive  with  a  warm  center  (higher  than  0.4><109 
J/m2)  located  in  the  southeast  part.  The  tropical  North  Atlantic  Ocean  (0°  -  25°N)  is  generally 
negative  with  a  cold  center  (lower  than  -0.6  x10  J/m")  near  Brazilian  coast.  The  mid-latitude 
North  Atlantic  Ocean  (25°  -  50°N)  is  generally  positive  with  a  warm  center  (higher  than  1.3  x109 
J/m")  near  Newfoundland,  Canada.  The  high-latitude  North  Atlantic  Ocean  (north  of  about 
50°N)  is  generally  negative  with  a  cold  center  (lower  than  -0.5  x109  J/m2)  near  Greenland.  The 
second  principal  component, PC2(t  )  ,  fluctuates  between  -1  to  0.72  during  1990-2009,  with  no 

evident  trend.  Since  the  EOF-2  mode  only  contains  8.83%  of  the  total  variability,  it  will  not  be 
discussed  here. 

6.  Conclusions 

The  major  goal  of  this  study  was  to  identify  the  interannual  variability  in  the  upper  OHC 
anomaly  and  its  relationship  to  climate  change.  A  3D  global  gridded  ( T ,  S)  dataset  has  been 
established  with  monthly  increment,  1  °x  1 0  horizontal  resolution,  and  33  vertical  layers  as  the 
World  Ocean  Atlas  (http://www.nodc.noaa.gov/OC5/WOAQ5/pr  woa05.html)  from  the  global 
temperature  and  salinity  profile  program  (GTSPP)  data  using  the  OSD  method.  The  upper  OHC 
for  the  global  oceans  was  calculated  from  the  3D  global  gridded  dataset  with  the  upper  layer 
thickness  as  300  m  for  the  Pacific  and  Atlantic  Oceans  and  as  150  m  for  the  Indian  Ocean.  The 
composite  and  EOF  analyses  were  conducted  on  the  OHC  data.  The  following  results  were 
obtained: 
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(1)  The  total-time  mean  OHC  field  (H  )  was  established  with  OHC  less  than  1.2  xlO9 
J/m2  north  of  40°N  and  south  of  40°S.  The  high  H(r.)  areas  are  associated  with  the  subtropical 

gyres  in  the  Atlantic  and  Pacific  Oceans  with  OHC  ranging  from  2.8  to  3.0  xlO10  J/m2.  The  OHC 
is  less  outside  than  inside  of  the  subtropical  gyres  in  the  same  latitudinal  zone,  such  as  in  the 
eastern  Pacific  (around  2.0xl010  J/m2). 

(2)  The  composite  analysis  indicates  the  following  features  of  the  mean  seasonal 
variability  of  OHC  (i.e.,  annual  mean  OHC  relative  to  the  total-time  mean  OHC).  The  maximum 

seasonal  variability  is  an  order  of  magnitude  smaller  than  the  total-time  mean  H  ( r ) .  Positive 

variability  (summer)  shows  eddy-like  structure  and  negative  variability  (winter)  is  more  spatially 
uniform. 

(3)  The  EOF  analysis  was  conducted  to  analyze  the  OHC  anomaly  data  (removal  of 
total-time  mean  and  mean  seasonal  variability)  for  the  Pacific,  Indian,  and  Atlantic  Oceans  to 
identify  low-frequency  variability  of  the  OHC  anomaly.  The  patterns  of  the  first  two  EOF 
modes  are  different  for  the  three  oceans:  (a)  two  zonal  modes  in  the  tropical  Pacific  with  the 
EOF-1  (51.46%  variance)  representing  the  canonical  El  Nino/La  Nina  events,  and  the  EOF-2 
(13.71%  variance)  representing  the  pseudo  El  Nino/La  Nina  events,  (b)  one  zonal  mode  in  the 
tropical  Indian  Ocean  with  the  EOF-1  (24.27%  variance)  representing  the  basin-scale 
cooling/wanning  events,  and  the  EOF-2  (20.94%  variance)  representing  the  Indian  Ocean 
Dipole  mode,  and  (c)  no  zonal  model  in  the  tropical  Atlantic  with  the  EOF-1  (49.26%  variance) 
representing  the  basin-scale  cooling/warming  events,  and  the  EOF-2  (8.83%  variance, 
insignificant  mode). 

(4)  In  the  Pacific  Ocean,  the  EOF-1  mode  is  strongly  related  to  the  Southern 
Oscillation  with  a  maximum  correlation  at  1 -month  lead  of  the  EOF-1;  and  the  EOF-2  mode 
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(pseudo  El  Nino)  is  related  weakly  to  the  Southern  Oscillation  with  a  maximum  correlation  at 
no-lag  to  3-month  lag  of  the  EOF-2.  In  the  Indian  Ocean,  the  EOF-1  mode  (basin-scale 
cooling/wanning  mode)  has  no  correlation  with  (SOI,  DMI).  The  EOF-2  mode  has  evident 
negative  correlation  to  SOI  with  a  minimum  value  of  -0.57  for  1 -month  lead  of  PC2 ,  and  evident 

positive  correlation  to  DMI  with  a  maximum  value  around  0.52  for  2-month  advance  of  PC2. 
This  clearly  shows  the  co-existence  of  the  positive  phase  of  the  EOF-2  mode,  the  negative  phase 
of  SOI,  and  the  positive  IOD  events.  This  implies  that  the  upper  ocean  may  take  a  key  role  in 
climate  variability. 

Appendix  A.  Optimal  Spectral  Decomposition 

A  new  data  analysis/assimilation  scheme,  the  optimal  spectral  decomposition  (OSD),  has 
been  recently  developed  to  analyze  fields  from  noisy  and  sparse  oceanographic  data.  Let  (x,  z)  be 
horizontal  and  vertical  coordinates  and  t  be  time.  A  physical  variable  c(x,z,t )  at  depth  zk  is 
decomposed  using  the  generalized  Fourier  series  (Chu  1999;  Chu  et  al.  2003a,  b,  2005a,  b; 
Ivanov  and  Chu  2008) 

M 

c(x,zi,0  =  c0(zJt,0  +  2]am(z/fc,t)  Om(x,z/t),  xeR(zk)  (Al) 

m= 1 

where  Co  is  the  horizontal  mean  of  the  variable  c.  M  is  the  truncated  mode  number,  Om(x,zt) 
and  am(zk,t )  are  the  orthogonal  basis  functions  (or  called  modes)  and  the  spectral  coefficients, 
respectively;  R(zk)  is  the  area  bounded  by  the  lateral  boundary  r ( zk )  at  depth  zk.  The  basis 
functions  { Om(x,zA  ) }  are  eigen- functions  of  the  horizontal  Laplace  operator 
( y-  ee  8"  Idx  +d  /  dy~  )  with  the  basin  geometry  and  certain  physical  boundary  conditions.  For 
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temperature  and  salinity,  the  homogeneous  Neumann  boundary  condition  is  taken  at  the  solid 
boundary  r(z)  (i.e.,  no  heat  and  salt  fluxes), 


=  -U„ 


-  =  0,  m  =  l,2,...,M, 


(A2) 


where  n  is  the  unit  vector  nonnal  to  T(z)  .  The  basis  functions  {®m  }  are  independent  of  the 

data  and  therefore  available  prior  to  the  data  analysis.  The  OSD  method  has  two  important 
procedures:  optimal  mode  truncation  and  determination  of  spectral  coefficients  {am}.  The 
optimal  mode  truncation  (M)  is  detennined  using  Vapnik  (1982)  variational  method.  Application 
of  the  generalized  Fourier  series  expansion  (1)  to  the  observational  points  with  P  as  the  total 
number  of  observations  leads  to  an  algebraic  equation 

Aa  =  QY .  (A3) 

where  a  =  (ai,  a2,  ...,  aM),  is  the  state  vector  (M-dimensional)  ;  A  is  a  P  x  M  matrix;  Q  is  a  P 
x  P  square  matrix  (P  >  M);  Y  is  a  /-’-dimensional  observation  vector,  consisting  of  a  signal  Y 
and  a  noise  Y\  Due  to  high  level  of  noise  contained  in  the  observations,  the  algebraic  equation 
(3)  is  ill-posed  and  needs  to  be  solved  by  a  rotation  matrix  regularization  method  (Chu  et  al. 
2004)  that  provides:  (a)  stability  (robustness)  even  for  data  with  high  noise,  and  (b)  the  ability  to 
filter  out  errors  with  a-priori  unknown  statistics.  Interested  readers  are  referred  to  Chu  et  al. 
(2003  a,  b;  2004). 
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Table  1.  Variances  of  the  first  six  leading  EOFs  of  H300  for  the  Pacific  Ocean. 


EOF 

Variance 

Cumulative  Variance 

1 

0.5146 

0.5146 

2 

0.1371 

0.6517 

3 

0.0644 

0.7161 

4 

0.0473 

0.7634 

5 

0.0375 

0.8009 

6 

0.0307 

0.8316 

Table  2.  Variances  of  the  first  six  leading  EOFs  of  H150  for  the  Indian  Ocean. 


EOF 

Variance 

Cumulative 

Variance 

1 

0.2427 

0.2427 

2 

0.2094 

0.4521 

3 

0.1172 

0.5648 

4 

0.1107 

0.6755 

5 

0.0548 

0.7303 

6 

0.0475 

0.7778 

Table  3.  Variances  of  the  first  six  leading  EOFs  of  H300  for  the  Atlantic  Ocean. 


EOF 

Variance 

Cumulative 

Variance 

1 

0.4926 

0.4926 

2 

0.0883 

0.5809 

3 

0.0561 

0.6370 

4 

0.0514 

0.6884 

5 

0.0401 

0.7285 

6 

0.0362 

0.7647 
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DMI 


Fig.  1.  Time  series  of  (a)  SOI  and  (b)  DMI. 
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Fig.  2.  Temporal  distribution  of  GTSPP  stations  during  1990-2009  for  this  study. 
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Fig.  4.  Total-time  mean  of  global  upper  ocean  heat  content  during  1990-2009:  (a)  H300(r.)  ,  and 
(b)  S150(r). 
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Fig.  5.  Mean  seasonal  variability  of  upper  layer  heat  content  (unit:  J/m  )  (a)  H300(r/,t/)  and  (b) 
Hi^Mz) ,  with  upper  panels  for  January  and  lower  panels  for  July. 


30 


Latitude 


32 


PC2 

A 

1  - 


Fig.  8.  Trajectory  in  the  phase  space  (PCi,  PC2)  of  the  OHC  anomaly  H300  in  the  Pacific  Ocean 

from  January  1990  to  December  2009  (a  star  enclosed  by  a  circle).  Here,  the  number  enclosed  by 
a  circle  denotes  January  of  the  year  (90  means  1990,  ...);  the  arrows  show  the  time 
advancement;  the  red  color  refers  to  negative  values  of  SOI  and  blue  color  refers  to  positive 
values  of  SOI.  It  is  noted  that  both  positive  (blue)  and  negative  (red)  SOI  are  found  in  all  four 
quadrants. 
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Fig.  9.  Lag  correlation  coefficients  between  (a)  SOI  (b)  DMI  with  PCi  (upper  panels)  and  PC2 
(lower  panels)  of  the  OFIC  anomaly  H300  in  the  Pacific  Ocean  with  positive  time  denoting  lag  of 
SOI  and  DMI.  It  shows  evident  correlation  between  SOI  and  PCi  with  1-month  advance  of  PCi, 
weak  correlation  between  SOI  and  PC2  with  no-lag  to  3-month  lag  of/JC2(/~;)  to  SOI,  and  no 

correlation  between  DMI  and  (PCi,  PC2).  Here,  the  two  dashed  lines  show  the  critical  values  for 
significant  lag-correlation  with  the  significant  level  of  0.05. 
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Fig.  10.  (a)  EOF-1,  and  (b)  EOF-2  of  the  OHC  anomaly  H150  (unit:  10s  J/m2)  in  the  Indian 

Ocean.  It  is  noted  that  the  EOF-1  mode  represents  basin-scale  wanning/cooling;  and  the  EOF-2 
mode  represents  the  dipole  mode. 
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Fig.  13.  Trajectory  in  the  phase  space  (PCi,  PC2)  of  the  OHC  anomaly  H150  in  the  Indian  Ocean 

from  January  1990  to  December  2009  (a  star  enclosed  by  a  circle).  Here,  the  number  enclosed  by 
a  circle  denotes  January  of  the  year  (90  means  1990,  ...);  the  arrows  show  the  time 
advancement;  the  red  color  refers  to  negative  values  of  DMI  and  blue  color  refers  to  positive 
values  of  DMI.  It  is  noted  that  both  positive  (blue)  and  negative  (red)  DMI  are  found  in  all  four 
quadrants,  but  more  in  the  third  and  fourth  quadrants. 
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Fig.  14.  Lag  correlation  coefficients  between  (a)  SOI  (b)  DMI  with  PCi  (upper  panels)  and  PC2 
(lower  panels)  of  the  OHC  anomaly  H150  in  the  Indian  Ocean  with  positive  time  denoting  lag  of 

SOI  and  DMI.  It  shows  no  correlation  between  (SOI,  DMI)  and  PCi,  evident  negative  correlation 
between  SOI  and  PC2  with  1 -month  advance  of  PC2,  and  evident  positive  correlation  between 
DMI  and  PC2  with  2-month  advance  of  PC2.  Here,  the  two  dashed  lines  show  the  critical  values 
for  significant  lag-correlation  with  the  significant  level  of  0.05. 
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(a)  EOF  Mode  1  (49.3%)  (xleP  J/m2)  (b)  EOF  Mode  2  (8.83%)  (xlC?  J/m2) 
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Fig.  15.  (a)  EOF-1,  and  (b)  EOF-2  of  the  OHC  anomaly  H300  (unit:  108  J/m2)  in  the  Atlantic 
Ocean.  It  is  noted  that  there  is  no  zonal  dipole  pattern  in  the  EOF-1  and  EOF-2  modes. 


